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THE KRIGIFIER: A PROCEDURE FOR GENERATING PSEUDORANDOM 
NONLINEAR OBJECTIVE FUNCTIONS FOR COMPUTATIONAL 

EXPERIMENTATION 

Michael W. Trosslt 1 

Abstract. Comprehensive computational experiments to assess the performance of algorithms for numer- 
ical optimization require (among other things) a practical procedure for generating pseudorandom nonlinear 
objective functions. We propose a procedure that is based on the convenient fiction that objective functions 
are realizations of stochastic processes. This report details the calculat ions necessary to implement our pro- 
cedure for the case of certain stationary Gaussian processes and presents a specific implementation in the 
statistical programming language S-PLUS. 
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1. Introduction. It is widely accepted that the performance of algorithms for numerical optimization 
should be established in fact as well as in theory. Factual evidence includes the anecdotal experiences of 
users, but it should also include (as do other empirical sciences) the results of carefully designed experiments. 
Unfortunately, it is not at all clear how to design meaningful computational experiments for numerical 
optimization. This report attempts to address that concern. 

Individuals who study numerical optimization often recommend specific algorithms for specific applica- 
tions. Typically, such recommendations are based partly on theory, partly on knowledge that the recom- 
mended algorithm has performed well on other, related applications. The latter rationale implicitly assumes 
that the relevant population of applications has been sufficiently well sampled to warrant making predictions 
about the new application in question. Is this usually the case? 

Computational experiments designed to assess the performance of algorit hms for numerical optimization 
have traditionally used a small number of canonical test problems. Most of these problems were created 
or discovered because they exhibit some special sort of pathology. Thus, t he fundamental premise of most 
computational experiments for numerical optimization is the following: the performance of an algorithm in 
typical situations can be inferred from its performance in certain pathological situations. Sadly, this premise 
seems dubious at best. 

Consider, for example, the simplex algorithm(s) for linear programming. In theory, the computational 
complexity of these algorithms is exponential; in practice, they invariably perform as if their complexity was 
polynomial. This discrepancy between worst-case and average-case performance had led some researchers 
to initiate theoretical studies of (j'ptctrd simplex performance on some simple populations of randomly 
generated linear programs. Although realistic distributions of linear programs undoubtedly render theoretical 
investigations intractable, one might still study empirical simplex performance on such populations. 
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As difficult as it may be to randomly generate plausible linear programs, it seems far more difficult to 
randomly generate plausible nonlinear programs. This report addresses one facet of the problem of generating 
random nonlinear programs, viz. the problem of generating random nonlinear objective functions. 

2. Basic Concepts. Originally developed by geostatisticians, kriging is a procedure for optimally 
interpolating a finite number of observed values of a realization of a specified stochastic process. (In case the 
stochastic process is an unspecified member of a specified parametric family of stochastic processes, kriging is 

preceded by estimation of the unspecified parameters.) The function / obtained by kriging values y\ , y T) 

observed at locations j’i,...,*,, is the expected value of the process, conditional on the process behaving 

as observed at x n . Thus. / can be regarded as a smoothed realization of the process and. <rt<ris 

pamhis , the degree of smoothing depends on how many values were observed: as more and more values are 
kriged, / looks more and more like an actual realization. 

As described in [1], the design and analysis of computer experiments is predicated on the fiction that the 
output from an expensive deterministic computer simulation resembles a realization of a stochastic process. 
VV e emphasize that this narrative is entirely fictional, convenient because it suggests plausible designs and 
analyses; nevertheless, simulations of complex physical phenomena often produce approximat ion, rounding, 
and truncation errors that contaminate the idealized output. Such deterministic noise can indeed resemble a 
realization of a stochastic process, so that it seems perfectly reasonable to synthesize inexpensive functions 
that approximate expensive simulation outputs by generating realizations of a stochastic process and adding 
each realization to a prescribed trend. 

( onceptually. the krigifier comprises the following steps: 

1. The user specifies an underlying trend, e.g. a quadratic function. 

2. The user specifies a stochastic process, e.g. a stationary Gaussian process. 

•h A finite number of points, are chosen at which the stochastic process will be observed. 

These points can be specified by the user or randomly generated by the krigifier. 

1. I he krigifier generates t/i, . . . , t/ n < the values of the stochastic process at i*i, .... ,r n . 

5. The krigifier interpolates y\, y n to obtain a noise term. 

(>. The trend and noise terms are added to produce an objective function. 

The next section describes each of these steps. 

3. Computational Details. This section describes precisely how the krigifier generates a pseudorandom 
nonlinear objective function: an implementation in t he statistical programming language S-PLIJS is provided 
in t he following section. 

1. Trend. A funct ion trend (a* ) is specified by the user. This might be a constant , e.g. trend (.r) = 0, but 
it seems more sensible to induce some underlying structure appropriate for nonlinear optimization, perhaps 
by specify ng a convex quadratic function. 



2. Stochastic Process, A stochastic process is specified by the user. The process should lx* one from 
which it is reasonably easy to generate a realization. We have experimented with stationary (Gaussian 
processes with covariance functions of the form 

c(.s, f) = ir 2 r(.s\ / ), (l) 

where rr 2 is the constant variance of the process at any point (the process is iiomoschedastic) and the 
correlation function is of the form 

>•(*./)= 6(||.s - f||) (2) 

(the process is isotropic). Specifically, we have experimented with 

O(u) = exp(— (3) 

for o=2 and o = 1. The former choice results in smooth (f ,rx ) interpolations that seem better suited 
to generating “nice" objective functions; the latter choice results in jagged interpolations that seem better 
suited to simulating numerical noise. 

3. Selected Sites. The user must specify /?, the number of sites at which the stochastic process will be 

observed. The locations of the sites can be chosen by any method whatsoever. In our experiments, we have 
specified a rectangle and drawn from a uniform distribution on the rectangle. 

4. Observed Values. Assume that the specified stochastic process is of the form described above. Given 

.r i , let 

n = {r(.r,..r-j)\ 

be the n x n matrix of interpoint correlations. We need to generate y = (// \ y u ) f by sampling from an 

n-variate normal distribution with covariance matrix rr 2 R. To do so. we exploit the fact that 

Theorem 1 // c ~ A : ((), / ). Uun Az ~ Y(0, .4 4'). 

First, we generate v standard univariate normal random variates, z\ c T7 . Next, assuming that R is 

positive semidefinite. let R = l r DI' f be its singular value decomposition. Then, letting c = (~i c,,)'. 

Theorem 1 tells us to set 

y = a( D l/2 t ":. 

5. Interpolation. We interpolate bv kriging. Assuming that R is invertible, define r by the square system 
of linear equations Rr = y. The information needed to define the interpolating function is contained in 
jq. . . . . i\ and the correlation function ?*(■.*). 

(iiven .r. let 

'•(■'•I..'') 

>■(•*■) = • 

_ r(x„. j ) 

noise(j-) = v'r(.r). 




Then the interpolating function is 



6. Additive Noise. The proposed pseudorandom objective function is 

f(.r) = trend (.r) 4- noise(r). 

4. An Implementation in S-PLUS. This section exhibits S-PLUS functions that perform the calcula- 
tions detailed in Section 3. The function krigify, exhibited in Figure 1. produces the information needed to 
define the noise term in the pseudorandom objective function /. The function f .rand, exhibited in Figure 2, 
evaluates /, which is constructed by adding the noise to a user-specified quadratic trend. 

funct ion (a , b,n, alpha , theta, sigma2) 

{ 

# 

# [a,b] is a p-dimensional rectangle; 

# n is the number of sites to be selected; 

# alpha, theta, sigma2 are parameters of an isotropic 

# stationary Gaussian process. 

# 

tol <- le-007 
p <- length(a) 

X <- matrix (runif (n*p ,min=a ,max=b) , byrow=T, nrow=n, ncol=p) 

R <- matrix(0, nrow=n, ncol=n) 
for (i in 2:n) { 
for (j in 1 : i) { 

R[i,j] <- exp(-theta * (vecnorm(X[i ,] -X[j ,])) “alpha) 

> 

} 

R <- R + t(R) + diag(n) 

Rsvd <- svd(R) 

d <- Rsvd$d [Rsvd$d >= tol] 

k <- length(Rsvd$d) - length(d) 

d <- c ( 1/ sqrt (d) , rep(0 , t imes=k) ) 

y o matrix (rnorm(n , sd=sqrt (sigma2) ) , ncol=l) 

v <- Rsvd$u 7, *7, diag(d) %*% t(Rsvd$u) 7. *7. y 

return(list (X=X , v=v)) 

> 


Figure (: The S-PLUS function krigify. 

1° krigify. it is necessary to specify a p-dimensional rectangle [«,6], the number n of sites to be 
selected from [«, />], and the parameters (o. 6h a 2 ) of a stationary Gaussian process with an isotropic covariance 
function of the form specified by equations (I), (2), and (3). The sites x n are drawn from a uniform 

distribution on [u, b]. 

Once the output from krigify has been saved, e.g. bv the S-PLUS command 
> noise <- krigify (a, b,n .alpha, theta, sigma2) 
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then it can be supplied to the pseudorandom objective function f whenever a function value is requested. 
This is accomplished by the S-PLUS function f .rand, exhibited in Figure 2. The function f .rand has two 
arguments, the ,r at which a function value f(,v) is requested and a list of auxiliary parameter values that 
specify /. and it returns /(c). 

f unction(x , aux) 

# 

# x is a p-dim vector at which f is to be evaluated; 

# aux is a list: 

# aux$betaO is a scalar, 

# aux$betal is a pxl matrix, 

# aux$beta2 is a pxp matrix, and 

# aux$xO is a p-dim vector that specify the quadratic trend; 

# aux$alpha & aux$theta specify the correlation function; 

# aux$X is an nxp matrix and 

# aux$v is an nxl matrix outputted from krigify. 

# 

n <- nrow(X) 

r <- matrix (nrow=n , ncol=l) 
for (i in l:n) { 

r[i,l] <- exp(-aux$theta * ( vecnorm(X [i , 1] -x) ) ~aux$alpha) 

} 

x <- matrix(x-aux$xO , nrow=length(x) , ncol=l) 

q <- aux$betaO + t(aux$betai) x + t(x) aux$beta2 '/,**/, x 

ret urn (q + t(aux$v) r) 

> 


Figure 2: The S-PLTS function f .rand. 

To illustrate the use of krigify and f .rand, Figure W exhibits SPIT S code for generating a pseudoran- 
dom objective function f on [0, l] 2 , evaluating / on a grid, and displaying the resulting function values in a 
perspective plot. 

5. Conclusions. We invite the reader to experiment with the krigifier and discover parameter settings 
useful for his or her applications. In our view, no amount of discussion can substitute for personal experience. 
Nevertheless, the krigifier does exhibit certain characteristics that deserve mention. 

1. Suppose that trend(c) is constant so that /(.r) = c T noise(.r). By construction, noise(.r,) ~ yj and 
noise(.r) tends to intermediate values of y for x 0 {.r lt . . . , x tl }. Hence, t he global minimizer of j in 

[a, b] will either equal or be near the global minimizer of / in the finite set {.ri r,,}. Because it 

is generally quite difficult to construct functions with multiple local minimizers and know the location 
of the global minimizer, the krigifier would appear to be especially useful for constructing global 
optimization test functions. 
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> a <- c (0 , 0) 

> b <- 0(1,1) 

> noise <- krigif y (a , b , 200 , 1 , SO , 100) 

> betal <- matrix(0 , nrow=2 , ncol=l ) 

> beta2 <- 100*diag(2) 

> xO <- c (0 . 3 , 0 . 4) 

> aux <- list(beta0=50, betal=betal, beta2=beta2, x0=x0 , 

alpha=l, theta=S0, X=noise$X, v=noise$v) 

> x <- y <- (0 : 50)/50 

> z <- matrix (nrow=51 ,ncol=51) 

> for (i in 1:51) { 

+ for (j in 1:51) { 

+ z[i,j] <- f .rand(c(x[i] ,y[j]) , aux) 

+ > 

+ > 

> persp(x,y,z) 


Figure 3: Using krigify and f .rand. 

2. Our own experience with the krigifier suggesls that it is easier to construct functions with multiple 
local minirnizers in low-dimensional spaces than it is in high-dimensional spaces. To provide a heuristic 
explanation of this phenomenon, suppose that { , . . . , x n } C [«,/>] C form a rectangular grid. A 

local minimizer will be induced at the grid point if each of t lie random variates assigned to the 2/> 
grid points adjacent to j* exceeds the random variate assigned to x. Obviously, the probability of this 
occurring decreases as p increases. To t he extent that realizations of stochastic processes are indeed 
plausible models of objective functions, this insight, suggests that functions of many variables may be 
l< s* likely to have multiple local minirnizers than functions of few variables, an amusing reversal of the 
curse of dimensionality. 
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